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Abstract 



In this paper, we develop a numerical procedure for investigating the dy- 
namics of trapped Bose gases based on the ZGN theory. The dynamical 
equations used consist of a generalized Gross-Pitaevskii equation for the con- 
densate order parameter and a semiclassical kinetic equation for the thermal 
cloud. The former is solved using a fast Fourier transform split-operator 
technique while the Boltzmann equation is treated by means of A^-body sim- 
ulations. The two components are coupled by mean fields as well as collisional 
processes that transfer atoms between the two. This scheme has been applied 
to a model equilibration problem, dipole oscillations in isotropic traps and 
scissors modes in anisotropic traps. In the case of the latter, the frequen- 
cies and damping rates of the condensate mode have been extracted from the 
simulations for a wide range of temperatures. Good agreement with recent 
experiments has been found. 

PACS numbers: 03.75.Fi, 05.30.Jp, 67.40.Db 
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I. INTRODUCTION 



The achievement in 1995 of Bose-Einstein condensation (BEC) in trapped atomic gases 
initiated a period of intense experimental and theoretical research which continues 
at an ever-increasing pace. Due to the relative diluteness of the atomic vapours, these 
systems present an excellent opportunity to investigate the effects of Bose degeneracy in 
novel situations. For example, a rich array of dynamical behaviour can occur as a result of 
the simultaneous occurence of a condensate order parameter and thermal excitations. 

It is now well established that at very low temperatures the dynamics of Bose-condensed 
gases can be accurately described by the Gross-Pitaevskii (GP) equation for the condensate 
order parameter (see e.g. @]). However, despite considerable work over the past few years, 
it has proved far more difficult to extend the theory to include thermal excitations. One 
of the major challenges is to capture the dynamics of the non-condensed component in a 
self-consistent way, while still allowing for a feasible calculation. A possible approach to 
this problem was suggested by one of the authors (EZ) together with Griffin and Nikuni 
PJ^. In this formalism (which we will hereafter denote by ZGN) the thermal excitations are 
described semiclassically by a Boltzmann kinetic equation, which couples to the condensate 
through mean fields and collisions. 

In this paper we report on recent progress in solving the ZGN equations to study the 
dynamics of Bose-condensed gases in experimentally relevant regimes. We begin by briefly 
reviewing the formalism, and present analytical calculations of collision rates as a function 
of temperature. We then go on to discuss numerical calculations, where the kinetic equation 
is solved using A^-body simulations with collisions treated through Monte-Carlo sampling. 
Applications to thermal cloud quenching, dipole modes, and scissors modes ^ are discussed. 
We conclude by highlighting future research directions. 



II. ZGN THEORY 

The ZGN theory Q is based on the assumption of a Bose broken symmetry which leads 
to the macroscopic wavefunction $(r, t) = ('?/'(r,t)). Here, ip{r,t) is the Bose fleld operator 
and the angular brackets denote an average with respect to a nonequilibrium density matrix. 
The condensate density is then given by nc(r, t) = |$(r, The condensate wavefunction 
is found to satisfy the generalized GP equation 

^fl§-^^r, ^) = ( - I^V^ + V{r) + g[n,ir, t) + 2n(r, t)] - iR{r, t)^ $(r, t), (1) 

where V{r) is the external conflning potential, n(r, t) is the noncondensate (or thermal 
atom) density, and -R(r, t) is a non-Hermitian source term that will be deflned below. The 
interactions between atoms are taken to be local and leads to the mean fleld potential 
g[nc + 2n], where g = ATrfi^a/m, with a the s-wave scattering length. For T — 0, n and R 
are zero and one recovers the familiar nonlinear GP equation. The stationary ground state 
solution $0 gives the density rico = |$oP which for large A^ is well-approximated by the 
Thomas- Fermi result rico = (/^o ^y)/g @]- Here, /io is the ground state GP eigenvalue or 
chemical potential. 
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Small amplitude excitations of the T = GP equation can be obtained by writing $ = 
$0 + ^^ and linearizing with respect to the deviation S^. This procedure is equivalent to the 
Hartree-Fock-Bogoliubov (HFB) approximation and generates the Bogoliubov quasiparticle 
excitations. At low energies, these excitations can be interpreted as collective excitations of 
the condensate, while at high energies they correspond to single-particle thermal excitations. 
However, at nonzero temperatures, one must take into account the occupation of excited 
states and the appearance of a noncondensate density n which interacts with the condensate. 
This problem was addressed by Hutchinson et al. (hereafter referred to as HZG) and Dodd 
et al. [0 who calculated tt-co and no self-consistently in the HFB-Popov approximation. The 
quasiparticle excitations calculated within this scheme become temperature dependent and 
at first sight might be associated with collective excitations of the trapped gas. However 
these excitations in fact correspond to collective excitations of the condensate in the presence 
of a static thermal cloud. In other words, the dynamics of the latter are not included 
consistently. This was revealed in calculations of the center of mass dipole mode whose 
frequency showed a deviation from the trap frequency cuq, which is the correct result as 
dictated by the generalized Kohn theorem |T^. This was the motivation for developing 
a theory which consistently included the dynamics of the thermal cloud with that of the 
condensate. 

The ZGN theory provides such a description. It is based on the following approx- 
imations: (i) the Popov approximation which neglects the anomalous average m(r, t) = 
{ip{r,t)ip{r,t)); here, ipir.t) = ip{r,t) — {ip{r,t)) is the noncondensate field operator, (ii) 
the use of Hartree-Fock thermal excitations with energy ep{r,t) = p'^/2m + U{r,t), where 
U = V + 2gn is the mean field potential acting on the thermal atoms, (iii) a semiclassical 
approximation for the thermal atoms which are represented by a phase-space distribution 
/(p, r, t), and (iv) the inclusion of binary collisions among thermal atoms as well as with 
the condensate. 

Within these approximations, one obtains the semiclassical Boltzmann kinetic equation 
for the thermal cloud 

^ + P . V/ - Vf/ ■ Vp/ = Cu[f] + C22[/]. (2) 
ot m 

The collision integrals appearing on the right hand side of this equation are given by 

Ci2[f] = ^^ j dp2 j c^Ps j dYii5{m\c + P2-Vz-Vi)K^c + e2-e^-ei) 

x[5(p - P2) - (5(p - P3) - 5(p - P4)] 

x[(l+/2)/3/4-/2(l + /3)(l + /4)], (3) 

C22[/] = ^f^^2 j j ^P3 j dpiSip + P2 - P3 - P4.)5{e + €2-63- £4) 

x[(l + /)(l + /2)/3/4-//2(l + /3)(l + /4)] . (4) 

Here, a = Svra^ is the total cross-section. The C22 integral involves the scattering of two 
atoms from initial to final thermal states. On the other hand, C12 describes collisions in 
which a condensate atom is involved in either the incoming or outgoing channels. These 
collisions have the effect of transferring atoms from or to the condensate and lead to a change 
in the number of noncondensate atoms, N. Specifically, we have 
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(5) 



dt J {2T^hY 

These collisions also define the source term 

appearing in Eq. ([l|). The condensate number A^^ of course satisfies dN^/dt = —dN/dt so 
that the total number of atoms, N = + N , is conserved. The local rates at which atoms 
enter and leave the condensate are given by 

= ^^3^2 J j ^P3 j 5(mVc + P2 - P3 - Pi)5iec + €2-63- £4) 

x(l + /2)/3/4. (7) 

x/2(l + /3)(l + /4). (8) 

In equilibrium, Vc = and Ec = fJ^o, the GP eigenvalue. One can then show that these two 
rates are equal if the distribution function takes its equilibrium Bose-Einstein form 

where (3 = l/Zc^T. As would be expected, the condensate and noncondensate must have 
the same chemical potential if the net transfer rate is to be zero. However, for the gen- 
eral nonequilibrium situation the rates are different and result in a time-dependent 
amplitude of the condensate wavefunction. This is revealed most clearly by making use of 
phase and amplitude variables, $(r,t) = a/^c(i", t)e'^*^'"'*^. The generalized GP equation is 
then equivalent to the continuity equation 

^ + V ■ (n.vj = rr^ - r°^* , (10) 

where Vc = kVO/m, and the Euler equation 

+ V (l;mv^] = -V/ic (11) 



dt \2 

with 



/ic = - +V + gn, + 2gfi. (12) 



It is instructive to examine r°2* for an equilibrium state of a typical trapped gas. r°2* 
represents the average number of thermal atoms colliding with the condensate per unit 
time per unit volume. Thus, dividing r°2* by the local noncondensate density gives the 
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local equilibrium collision rate per atom, l/'rj'g, for this kind of collision [|TT|. In Fig. 1 we 



plot '\-/t^2 fo'^ R-t) atoms in an isotropic trap with frequency z/q = 200 Hz for various 
temperatures below the ideal gas critical temperature T° ~ Q.^Ahuj^N^/^ /kB- The collision 
rate has a maximum at the edge of the condensate where the local fugacity z{r) = e^^'^°~^°^ 
approaches unity. Since 1/t^2 is proportional to rico, it falls off rapidly beyond this point. It 
is also interesting to note that for T approaching T°, ujqTi2 < 1, indicating that a thermal 
atom in this region of the condensate has a high probability of making a C12 collision on the 
timescale of the trap period. However, the number of thermal atoms which overlap with the 
condensate is relatively small and it is therefore more meaningful to define an average C12 
collision rate per thermal atom by 

"<^'rM- (13) 



ri2(T) N{T) J t\ 



This quantity is plotted in Fig. 3 which shows that the thermal cloud is in the transition 
regime ujoT^2 — 1 rather than the truly hydrodynamic regime co'oTi2 -C 1 as the local collision 
rate in Fig. 1 suggests. Nevertheless, these results show that C12 collisions play an important 
role in the dynamics of the thermal cloud in the region of the condensate. Of course, the 
importance of these collisions increases with increasing N and one can expect to be in the 
hydrodynamic regime for > 5 x 10^ for this particular trap configuration. 

One can also define an equilibrium collision rate between atoms in the thermal cloud. It 
is this rate that one would conventionally consider when deciding between coUisionless and 
hydrodynamic behaviour. Accordingly, we define a C22 collision rate through the relation 

-0 _ f dp, 



A J (Mr 

= /^''^ / ^''^ / ^^^^^1^2(1 + /3°)(1 + f!) ■ (14) 

Here, CI2* denotes the term in Eq. (^) in which particles in states '1' and '2' scatter into 
the final states '3' and '4'. In equilibrium, the C22* and C22 terms are equal. In Eq. (|T4|), 
Vj, = Vi — V2 is the relative velocity of the incoming particles and the solid angle fl defines 
the orientation of the final relative velocity vector, = V3 — V4. Momentum and energy 
conservation ensures that V1 + V2 = V3 + V4 and Vr = v'^. In the classical limit {h — >• 0), I/T22 
reduces to (1/t22)c/ = V^o'Vthno, where Vth = {SkT/nmy^'^. In Fig. 2 we plot I/T22 for the 
same trap parameters used to generate ^/t^2- -^^^ T <T^, we again see Bose enhancement 
of the C22 collision rate at the edge of the condensate, with ujqT22 ^ 1 being reached in 
this region at the higher temperatures. In Fig. 3 we present the average C22 collision rate 
1/070X22. Below T°, this rate is comparable to the average Cu collision rate, indicating that 
the probability of a thermal atom colliding with another thermal atom is similar to the 
probability of colliding with the condensate. However, unlike l/ooofu which involves the 
condensate density rico, l/u;or22 continues to increase up to T°, at which point it reaches 
a value of a;of22 — 1- This limiting value increases with increasing N, indicating that the 
hydrodynamic regime is easily reached near T°. Above T°, the effect of Bose enhancement 
rapidly diminishes and the average rate approaches its classical limiting value. 
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III. NUMERICAL METHODS 



Eqs. (D) and (|) can in principle be used to describe the dynamics of a trapped Bose gas 
for arbitrary conditions of temperature, number of atoms and trap geometry. The equations 
have been solved previously for various limiting situations. For example, Eq. (|l]) was solved 
by HZG 1^ for R = and with the assumption of a static thermal cloud. As stated earlier, 
this neglects the perturbation of the noncondensate that arises from the time-dependent 
mean-field of the oscillating condensate. As a result, it neglects those modes, such as the 
out-of-phase dipole mode, which involve the collective motion of the thermal atoms. Eqs. (|I]) 
and (H) have also been used to study collective modes in the collisionless limit using what 
is effectively a moment method |T^. In this approach, an ansatz for the form of the GP 
wavefunction and distribution function of the thermal cloud is made, and a finite number of 
moments of the GP and kinetic equations are retained to obtain a closed set of equations for 
the modes of interest. However, due to the truncation made, the effect of Landau damping on 
the modes is not included. Furthermore, the ansatz used for the form of the nonequilibrium 
distribution function does not treat the mean-field interactions between the two components 
accurately. In the opposite limit of high collision rates, we previously developed various 
approximations which address the dynamics in the hydrodynamic regime 0. By assuming 
C22 collisions to be rapid we developed a description in which the thermal cloud is assumed 
to be in local thermodynamic equilibrium, but not necessarily in local equilibrium with 
the condensate The assumption of complete local equilibrium results in a Landau- 
Khalatnikov two-fluid theory 



Clearly it would be desirable to obtain accurate solutions of our equations which avoid 
approximations beyond those used to generate the equations themselves. In this way, the 
limitations of the approximations underlying the theory can be critically examined. For this 
purpose we have developed a numerical procedure within which the accuracy of the solutions 
can be studied and evaluated. We outline our numerical methods in the following. 

The generalized GP equation can be solved by a split-operator technique. For the pur- 
poses of the following discussion, we write Eq. (|l]) in the operator form 

(9$ 

^h— = {T + V)<!>, (15) 

where T is the kinetic energy operator and K is a complex, time-dependent potential. The 
wavefunction is advanced in time using the split-operator method: 

Since the V operator is local in position space the effect of e~*^'^*/^'^ is most easily treated 
in the coordinate representation. The kinetic energy operator on the other hand is local in 
momentum space and it is therefore useful to Fourier transform the intermediate state when 
determining the effect of e~*-^^*/^. A fast fourier transform (FFT) routine is used for this 
purpose. The inverse transform is then applied in order to complete the evaluation of the 
evolution operator. Since V is time- dependent, a higher order approximant to the evolution 
operator is in fact used to increase the accuracy and stability of the GP equation solution. 

The use of the semiclassical Boltzmann equation for the thermal atoms corresponds to 
a Hamiltonian description of the particle dynamics. In our case the Hamiltonian is simply 
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Hci = p^/2m + U{r,t), where U{r,t) = V{r) + 2gn{v,t) is the time-dependent mean-field 
potential. Thus the trajectory of a given atom between collisions can be determined using 
Newton's equations of motion. The evolution of the distribution function /(p, r, t) is then 
determined by considering the dynamics of Nf representative test particles |T5Hl7| which in 
effect corresponds to the distribution function /r(p, r, t) = Xlil^i ^(^ ~ i"j)'^(P ~ Pi)- the 
absence of collisions, the dynamics of this swarm of test particles is equivalent to solving the 
collisionless Boltzmann equation. 

In reality, the deterministic evolution of a particular phase space trajectory is interrupted 
stochastically by collisions with other atoms or with the condensate. In other words, in a 
time step At, a given thermal atom has some probability of suffering a collision, and the 
numerical algorithm used to generate these probabilities must yield collision rates that are 
consistent with the collision integrals in the Boltzmann equation. To simulate the C22 
collisions |T8|, we bin the test particles into position space volume elements A^r. Pairs of 
atoms (ij) chosen at random from a given volume element define collision partners, and the 
probability of their suffering a collision in a time interval At is given by 

P,,=na|v,-v,.| J ^(l + /3)(l + /4)At. (17) 

The angular integral is an average over all possible orientations of the final relative velocity 
and can be replaced by a single scattering event by choosing a random solid angle Qr 
uniformly distributed on a unit sphere. Thus the probability used in the simulation is 

Pi, = ha\v, - v,|(l + /3^-)(l + /r«)At . (18) 

The final velocities V3 and V4 are determined by Vj, \j and through momentum and 
energy conservation. To test whether a collision actually occurs, a random number Rij 
distributed uniformly between and 1 is chosen and compared to Pij. If Rij is less than 
Pij, the collision is accepted. This process is repeated for each pair (ij) in every real-space 
volume element. We have checked that the procedure correctly reproduces the equilibrium 
collision rates as given by Eq. (0). 

The C12 collisions are treated in a similar fashion. In this case, both 'in' and 'out' 
collisions, as illustrated schematically in Fig. 4, must be considered separately since one of 
the incoming or outgoing particles is from the condensate. The overall rate of 'out' collisions 
is given by Eq. (|]), in which a particle with velocity V2 collides with the condensate to 
produce two outgoing thermal atoms with velocities V3 and V4. This process is similar to 
a C22 collision in that the final relative velocity = V3 — V4 can have arbitrary direction 
with a magnitude of v[. = a/|vc — V2P — Agric/m. This value is a consequence of energy 
conservation and accounts for the fact that the mean-field energy of a thermal atom differs 
from that of a condensate atom by gric- Choosing one of the test particles, i, from the 
volume element A^r, the probability of it suffering an 'out' collision is given by 

Pr' = n^avUl + /3^«)(1 + f2-)At . (19) 

The situation for 'in' collisions is somewhat more complex since an arbitrary pair of ther- 
mal atoms in A^r will be precluded by energy and momentum conservation from scattering 
into the condensate. To ensure that these collision events are sampled correctly, we rewrite 
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Eq. (0) by interchanging the '2' and '3' labels in the integral. The '2'-particle is then paired 
with the '4'-particle to produce the outgoing '3'-particle. The resulting integral leads to a 
significantly different form for the probability of this type of collision event. Without giving 
details, we find that {i = '2') 

Pr = ^^£^(1 + f;-)f:-At , (20) 

where is a random point in a plane of area in velocity space. The plane is defined by 
energy and momentum conservation, and once the velocity is chosen, the velocities V3 and 
V4 are determined. Pj™ is largely independent of Ay as long as the plane completely samples 
the occupied regions of phase space. We have checked that P°"* and P*" give identical 
collision rates for the equilibrium situation. The simulation of Cu collisions then proceeds 
as follows. For each particle i in A'^r, P°"* and P/"" are calculated. A random number Ri 
is selected and if Ri > P?" + P°"*, no collision occurs, while if Ri < Pl^, the 'in' collision 
is accepted and if P/" < Ri < P-"" + P°"*, the 'out' collision is accepted. This procedure 
ensures that any given particle only makes one kind of collision, if a collision occurs at all. 
We have again checked that this prescription gives correctly the equilibrium 'in' and 'out' 
collision rates. 

We finally address the question of determining the mean-field potential, 2gh, for A^^ 
point-like test particles. An approximation to n{r) is obtained by dividing space into cells 
A^r and binning the distribution to construct a density histogram. If A^r is sufficiently small 
and Nt sufficiently large, this procedure will generate a fairly smooth density distribution. 
Since our simulations are restricted practically to Nt < 10^ and a three-dimensional spatial 
grid that is relatively coarse, statistical fiuctuations in the number histogram are inevitable. 
To compensate for these fiuctuations we have adopted the following strategy. First, we use 



a cloud-in-cell method [jT9[ to define density weights at the 8 grid points bounding a cell 
according to the position of a particle within the cell. Having defined a density on the grid 
points of our mesh, we next convolve the density with a gaussian: 



7^3/2^3 



^^/g-|r-r'|Vr;2^(r') . (21) 



By choosing a width parameter rj somewhat larger than the spatial grid used in our calcu- 
lations, we achieve a further smoothening of the thermal atom density. The final mean-field 
potential of the thermal atoms is then defined as 2gnconv{^)- This result can also be viewed 
in terms of an interaction potential g{v — r') = (^f/vr^/^r/^) exp(— |r — r'p/77^) which has a 
finite range rj as opposed to the contact interaction g6{r — r') usually assumed. Such a 
procedure has a minor effect for density variations on length scales larger than rj but has 
the benefit of reducing fiuctuations stochastically generated on the scale of the spatial grid. 
We can of course check that the physical properties of interest are unaffected by the value of 
1] used in the calculations. The gaussian convolution has the further useful feature of being 
readily calculable with a FFT routine. 

Our simulations commence with the specification of the initial state of the system, namely 
the condensate wavefunction $(r,0) and the initial distribution of test particles. The par- 
ticular choice is dictated by the physical situation being treated and several examples are 
considered in the following section. Once the initial state is known, the condensate and 
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noncondensate mean-field potentials arc determined. The condensate wavefunction is prop- 
agated for a time At and the test particle positions and momenta are updated. At the 
end of this time step, the C22 and C12 collisions are considered in turn and their net effect 
on the thermal cloud distribution is determined. At the same time, the C12 coUisions are 
accumulated to construct the source term R{r,t) that is needed for the next cycle in the GP 
evolution. This procedure is then repeated for the duration of the simulation. Quantities 
of physical interest can be calculated from the condensate wavefunction and test particle 
distribution at each time step and are recorded for subsequent analysis and display. 



IV. APPLICATIONS 
A. Equilibration 

We begin by considering a relatively simple model problem which illustrates the equili- 
bration of an initial nonequilibrium state. We imagine a trapped gas in equilibrium at some 
temperature Tq which is chosen such that approximately half of the atoms are in the con- 
densate. To be specific, we consider 5 x 10^ ^^Rb atoms in an isotropic trap with harmonic 
frequency uq = 187 Hz. The scattering length is a = 5.82 nm. The equilibrium state is 
determined by solving self-consistently the following pair of equations: 

- ^^$o(r) + [V{r) + ^neo(r) + 2^no(r)]$o(r) = /xo$o(r) , (22) 

and 

no(r) = ^^3/2(^(r)) , (23) 

where A = (27r/i^/m/cTo)^/^ is the thermal de Broglie wavelength, zo{r) = exp(/3o(/io — 
Uo{r)) is the local equihbrium fugacity, Po — I/ZcbTq and Uo{r) — V{r) + 2g[nco{r) + 
no(r)]. For Tq = 200 nK, we find Nco = 2.58 x 10^. We next define a nonequilibrium state 
by first constructing a sample of test particles representative of the equilibirium density 
no(r). Each of these particles is then distributed in momentum according to the phase 
space distribution /(p, r) = (A/Ao)^[2;o(i')~^ exp(/5p^/2m) — where the temperature T 
is half the equilibrium value Tq. This initial nonequilibrium state has a total density equal 
to the equilibrium distribution, but the cloud is too 'cold' to sustain this distribution as 
a function of time. It begins to collapse in real space, with C12 collisions simultaneously 
transferring atoms from the thermal cloud to the condensate. In Fig. 5 we show the evolution 
of the condensate and noncondensate numbers as a function of time. One can see that the 
condensate number relaxes to a higher value corresponding to a final equilibrium temperature 
T' < Tq, on a timescale expected on the basis of the equilibrium C12 collision rates. The figure 
also shows the total number of atoms, = Nq + iV, as a function of time in the simulation. 
This number should of course be constant but the simulation does not conserve particle 
number exactly at each time step because of the different ways in which particle number 
is changed in the condensate and thermal cloud. The Monte-Carlo samphng of coUisions 
leads to integer- valued changes in the number of thermal test particles. On the other hand, 
the normalization of the condensate changes quasicontinuously at a rate determined by the 
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strength of the source term R{r, t) in the GP equation. The fact that in our simulations 
exhibits only small fluctuations about the initial value is reassuring and demonstrates that 
number conservation is preserved in an average sense. This is an obvious requirement of any 
method that hopes to address problems such as condensate growth ||2^ . 

In Fig. 6 we show the root-mean-squared radius a/ (r^) of both the condensate and 
thermal cloud. The thermal cloud initially begins to collapse as a result of the quench 
before bouncing back and oscillating at a frequency very close to the monopole frequency 
2ujq of a noninteracting thermal cloud. Oscillations in the condensate are in turn induced by 
the time-dependent mean-field potential as well as by the condensation of thermal atoms. 
These oscillations at first build up in time and then damp out as the entire system approaches 
overall equilibrium. The frequency of these oscillations is about 5% lower than the TF 
frequency for the condensate monopole mode, \f^ujQ This confirms that our simulations 
are providing a good description of the dynamical behaviour of both the condensate and 
thermal cloud. 



B. Dipole Oscillations 

Our second example is perhaps an even more stringent test of our numerical procedures. 
For harmonic confinement, a trapped gas must exhibit a centre of mass oscillation in which 
the total density of the gas oscillates rigidly without change in shape at the trap frequency 
a;o, i.e., n(r, t) = no(r — r}{t)), with r}(t) = r/o cos(to'ot). This result is a consequence of the 
generalized Kohn theorem |jTO| and must be satisfied in any consistent theory of the dynamics. 
The failure to satisfy this property was the signature that the HFBP mode calculations were 
deficient The ZGN theory, however, is formally consistent with the generalized Kohn 
theorem ^j. Our purpose here is to demonstrate that our numerical implementation of the 
theory satisfies this requirement as well. 

In these simulations, we retain the same trap parameters as given in the previous sub- 
section and the same number of Rb atoms, but at a higher temperature (Tq = 225 nK) 
where the difference between the condensate and noncondensate numbers is larger. The 
reason for wanting this difference to be appreciable will become apparent shortly. We begin 
by obtaining the equilibrium condensate wavefunction and an equilibrium distribution of 
test particles. To excite the dipole oscillations, we then displace each component a small 
distance Az along the 2;-axis relative to the trap, but in opposite directions. This is easily 
achieved by shifting the condensate wavefunction a certain number of grid points on the 
spatial mesh, and adding a constant displacement Az to the position of each test particle. 
The two components are then released and their subsequent motion is determined. 

In our first simulation, we turn off all coUisional processes. The two components then 
interact only as a result of the mean-field potentials. The mean displacement of the con- 
densate, {zc{t)) = J drznc{r,t) /Nc, and the thermal cloud, {z{t)) = J2f=i^ii'^)/^T, are then 
calculated as a function of time. In Fig. 7 we plot the center of mass position 

Zcmit) = ^ (iVc(^c(t)) + N{~z{t))^ , (24) 
and the relative displacement 

zit) = {z^it)) - {zit)) ■ (25) 
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We see that Zcmif) oscillates at the trap frequency, cuq, as expected, with no change in 
amplitude. This amplitude is given by Zcm{0) = Az{N — Nc)/N and is non-zero because 
N ^ N^in the equilibrium state. The relative displacement z{t) is also shown, and it is seen 
to decay as a function of time, although its behaviour at longer times is more complex. 

These results can be interpreted in terms of a two mode model 0: (i) a centre-of- 
mass mode for which {zc{t))i/ {z{t))i = 1 and (i) an out-of-phase dipole mode for which 
{zc{t))2l {z{t))2 = —N/Nc- Within this picture, the initial displacement of the two compo- 
nents excites a superposition of these two modes and the individual displacements are given 
by 



Az 



{z^(t)) = (^{N - N^)/N^ coscjit - {2N/N) cosc^a^ 

{z{t)) = (^{N - Nc)/N)^ coscuit + i2NjN) coscjat Az , (26) 

where Ui and uj2 are the two mode frequencies. In this model, Zcm{t) = Az{N — 
Nc)/N cosujit and z(t) = —Az{4N/N)cosuj2t. These expressions qualitatively describe 
the results of the full simulation, although the relative displacement is actually damped as a 
result of coupling to other internal degrees of freedom. As the condensate moves through the 
thermal cloud, the thermal atoms experience a time- dependent mean-field potential which 
can transfer energy to them. This excitation mechanism is the source of Landau damping 
of the mode. However, the interplay between the internal degrees of freedom and the out- 
of-phase collective variable is complex and is the reason why z{t) does not exhibit a simple 
damped-exponential behaviour. 

In our second dipole simulation, we include both the Cu and C22 collision processes. The 
qualitative behaviour is similar to our previous simulation in that the centre-of-mass mode 
continues to oscillate without damping at the trap frequency, Uq. This is the signature that 
the generalized Kohn theorem is being satisfied. In contrast to the centre-of-mass mode, the 
out-of-phase dipole mode is damped. Interestingly, the behaviour is more regular than in the 
absence of collisions. The damping is still due to the Landau mechanism, but as these mean- 
field excitations proceed, collisions are continually driving the gas towards local equilibrium. 
These collisions limit the extent to which the thermal cloud phase space distribution can 
be driven out of equilibrium, and as a result, the complex long-time behaviour seen in the 
absence of collisions is eliminated. 



C. Scissors Mode Oscillations 

As our final application, we consider the scissors mode oscillations that have recently 
been studied in anisotropic traps [^,^]. A preliminary report of our work has appeared 
and we shall here summarize some of the main findings. 

We consider an anisotropic harmonic trap of the form V{r) = miu'^pP^ + u;^z^)/2 which 
has axial symmetry about the 2;-axis, with ujz — \/9>ujp. In the experiments I^Sf, the scissors 
mode is excited by adiabatically rotating the trap potential through a small angle 9q and 
then suddenly rotating the trap through an angle —26q. To simulate this process, we first 
determine the equilibrium state of the system at some temperature T using the method out- 
lined earlier. We then rotate the particle coordinates and condensate wavefunction relative 
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to the trapping potential through an angle 26o about the y-axis. At this point the dynamical 
simulations commence. 

The relevant dynamical variables in this situation are the a;z-quadrupole moments 
Qc(t) = f dr xznc{r,t) and Q(t) = {N/NT)'^fjiXi(t)zi(t). If we assume that the two 
components rotate about the y-axis without distortion, these quadrupole moments would 
be given by Qa = 0(t)Q'^, where 9{t) is the instantaneous rotation angle and Q° = {x^ — z'^)'^ 
is the equilibrium quadrupole moment of the a-th component. We thus define the 'rotation 
angle' 9a(t) = Qait)/Qa as the variable representing the angular displacement of the two 
components in our simulations. The advantage of using the angular variables is that the 
range of angular displacements is similar for both components, and independent of temper- 
ature, in contrast to the quadrupole moments. 

In Figs. 8-10 we show the angular displacements of the two components as a function of 
time for a few temperatures. At the intermediate temperature of T = 185 nK, N^/N = 0.33, 
and we see the condensate exhibiting a damped oscillation at a single frequency, while the 
thermal cloud oscillation consists of a superposition of two frequencies. This is to be expected 
on the basis of the modes that are found at low temperatures for a pure condensate and for 
a thermal cloud above T° For the T = condensate, a straightforward analysis of the 



linearized TF equations of motion which follow from Eqs. (|TDp-(|T2|) allows one to determine 
the scissors mode frequency. Multiplying the continuity equation by xz and integrating over 
all space leads to 

^(x2;) = (xvcz + zvcx) , (27) 

where the velocity moments are defined as (xiVcj) = J dr XiVcj{r,t)nco{r). Taking these 
moments of the velocity equation gives 

^{xvcz + zVcx) = -{ujI + tol){xz) . (28) 

These equations define the condensate scissors mode with frequency ujsc = a/cu^ + cu^. Our 
T = simulations for N = 2 x 10^ yield a frequency which is 1.5% larger than the TF result 
due to finite number effects. 

A similar moment analysis can be carried out for the Boltzmann equation above T° p4 . 
One obtains the following moment equations 

^{XZ) = {xVz + ZVx) 

dt 

— (xw^ - zvx) = 2eujl{xz) 

^{xvz + zvx) = 2{vxVz) - 2ujI{xz) 

^{vxVz) = -ujI{xVz + zvx) - eujl{xvz - ZVx) - '^^-^-iL , (29) 
at T 

where e = {uj'^ — ujI)/{uj^+ujI) and Uq = {ojI+ojI) /2. In the collisionless limit (r — > oo), these 
equations yield two scissors modes with frequencies uj± = Uz ^ ujp. Our simulations above 
T° agree well with these predictions. As Figs. 8-10 indicate, these modes persist down to 
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quite low temperatures below T° due to the relatively weak mean-field interactions between 
the condensate and thermal cloud. 

We now discuss the temperature dependence of the scissors mode oscillations. Fig. 8 
shows results for T = 55 nK and NJN — 0.95. The condensate mode is only weakly 
damped with a frequency which is close to the T — limit. However, the dynamics of the 
thermal cloud is quite different from the behaviour seen at higher temperatures. At this 
low temperature, the dynamics of the thermal cloud is strongly affected by the condensate 
mean-field. The fact that the thermal cloud shows an oscillation at the condensate fre- 
quency is indicative of its being driven by the more massive condensate. With increasing 
temperature, this coupling weakens and the thermal cloud begins to display its own scissors 
mode frequencies (Fig. 9). At the same time, the damping of the condensate mode is seen 
to increase as a result of the increasing rate of Landau damping. Finally, at the highest 
temperature of 296 nK shown in Fig. 10, N^/N = 0.02, and the roles of condensate and 
thermal clould are reversed in that the condensate is strongly coupled to the more massive 
thermal cloud and exhibits a more complex oscillation pattern. 

The interaction between the condensate and thermal cloud in the course of their oscilla- 
tions is revealed by a density contour plot for the two components. This is shown in Fig. 11 
for a sequence of equally spaced time steps spanning approximately half a period of the con- 
densate oscillation. What is striking about this figure is the lack of any obvious distortion 
of the condensate. It appears very much to be a rigid body (but not in a rotational sense!) 
as it oscillates in the anisotropic trap. The thermal component on the other hand shows 
clear signs of its interaction with the repulsive mean-field of the condensate. The condensate 
tends to move the thermal cloud out of its way in the course of its motion, much like a boat 
moving through water. Clearly the thermal cloud is strongly perturbed in the vicinity of the 
condensate. It is here that excitation of the thermal atoms, which leads to Landau damping 
of the condensate oscillation, occurs most strongly. 

A comparison of our results with experiment is shown in Fig. 12. To obtain the fre- 
quencies and damping rates of the various scissors modes we make the assumption that the 
calculated time- dependent quadrupole moments are a superposition of three normal modes, 
one corresponding to the condensate scissors mode and the other two to the thermal cloud 
scissors modes. Accordingly, we fit the time-dependent data in Figs. 8-10 with three damped 
sinusoidal functions with adjustable amplitudes. The values of the frequency and damping 
rate of the condensate scissors mode extracted in this way are found to be in very good 
agreement with experiment for T < 0.8T°. Discrepancies occur above this temperature, 
but here the experimental error bars are quite large. As discussed above, the condensate 
is strongly coupled to the thermal cloud as T° is approached and it becomes difficult to 
obtain good fits to the condensate quadrupole moment using the three-mode analysis. This 
is probably also the source of the large error bars associated with the experimental fits. 

At low temperatures, the frequency of the condensate mode at first increases slightly 
with increasing temperature. This we believe is due to an increasing deviation from the TF 
frequency ^/ouj + uil as a result of the decreasing number of condensate atoms. At higher 
temperatures, the frequency begins to decrease as a result of the increasing size of the 
thermal cloud. As the density contours in Fig. 11 suggest, the condensate drags part of the 
thermal cloud along with it and one would expect the condensate to effectively have a larger 
inertia, and therefore a lower frequency. The damping of the condensate mode is essentially 
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due to Landau damping but Cu and C22 collisions are indirectly playing an important role. 
Without collisions, the damping rate is as much as 50% lower than the values shown in the 
figure. Collisions have the effect of equilibrating the thermal cloud which is being driven 
by the condensate, and it appears that this internal equilibration is enhancing the rate at 
which excitations of thermal atoms occur. 

The existence of a single condensate scissors mode is a conseqence of the irrotational 
nature of the condensate velocity field. In contrast, the velocity fields of each of the two 
thermal cloud modes have both irrotational and solenoidal components. This distinction 
becomes apparent in the moment of inertia of the system, /. Zambelli and Stringari 
have shown that the moment of inertia is given by the expression 

(30) 

/rigid j dux {uj) uj 

where /rigid is the moment inertia of the system treated as a rigid body. The quantity x"{^) 
is the imaginary part of the quadrupole response function of the system: 

xH = -^ dt{[Q{t),Qm)o, (31) 

where Q = YliLi ^i^i quadrupole moment operator. This function determines the 

linear response of the system to an arbitrary time-dependent perturbation that couples to 
Q. In particular, the sudden rotation of the trap through a small angle do induces the 
time-dependent quadrupole moment Q{t) determined in our simulations. It can be shown 
that 

x"iuj) = -uj^{Q'{iu)/x} , (32) 

where Q'{uj) is the Fourier transform of Q{t) — Q{oo) and A = m{uj'^ — ujI)6o- Our three- 
mode analysis of the induced quadrupole moment yields a x"{^) consisting of a sum of 
Lorentzian spectral densities for which the frequency moments required in Eq. ( pUD are 
undefined. However, the moment of inertia can be estimated by replacing the Lorentzian 
line shapes by delta functions having the same spectral weight. Using this procedure, we 
show the moment of inertia of the system as a function of temperature in Fig. 13. At high 
temperatures, / tends towards /rigid = m{x'^ + z'^)N, as expected when the large thermal 
cloud dominates the response. On the other hand, as T ^ 0, the moment of inertia takes 
on the irrotational value, / = /rigid, where e = (x^ — + 2;^)° parameterises the 

deformation of the condensate in equilibrium. 

Under the semiclassical approximation, Stringari p6[ showed (at least in the ideal gas 



case) that / at intermediate temperatures is simply the sum of the two limiting values, 
weighted by the number of atoms in each component 

I = fm{x^ + z^)cNc + m{x^ + z'^)nN . (33) 

This quantity can be readily calculated from the equilibrium densities and is plotted (as a 



ratio with /rigid) in Fig. |I3|. Comparison to the results from (^) show that both monoton- 
ically increase with rising temperature, clearly indicating the transition between superfiuid 
and normal fiuid behaviour. Small discrepancies between the two approaches may result 
from our evaluation of the frequency spectrum, or in generalising (0) to an interacting gas. 
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V. CONCLUSIONS 



We have used the ZGN theory to investigate the dynamics of trapped Bose gases at fi- 
nite temperatures. In this theory, the condensate is described by a generahzed GP equation 
while the thermal cloud is treated semiclassically in terms of a Boltzmann kinetic equation. 
We have developed a numerical procedure for solving these equations and have applied the 
method to a variety of dynamical problems including equilibration and collective excita- 
tions. We have demonstrated that the ZGN theory provides an effective method for dealing 
with such problems and is a useful tool in the interpretation of experiment. In particular, 
the method has been used to examine the scissors modes in anisotropic traps. The good 
agreement with experiment found for this entire class of modes is encouraging and suggests 
that other situations can also be analyzed successfully. We plan to contimic this work with 
a systematic study of the various collective modes that have been probed experimentally, as 
well as problems associated with vortex nucleation and dynamics. 
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FIGURES 
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FIG. 1. The equilibrium local C12 collision rate l/rj'a in units of the trap frequency, ujq, as 
a function of the radial distance in units of the harmonic oscillator length d = ^KjmujQ. The 
different curves are labelled by the temperature T. The trap is isotropic with harmonic frequency 

= 200 Hz and contains 10^ ^^Rb atoms. 
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FIG. 2. The equilibrium local C22 collision rate I/T22 plotted as 1/t]^2 1- The curves 

labelled 1 through 5 correspond to temperatures between 100 nK and 500 nK in steps of 100 nK. 
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FIG. 3. Average collision rates per atom in units of the trap frequency, wq, as a function of 
temperature. The solid squares are for C12 collisions and the open squares are for C22 collisions. 
The trap parameters are the same as in Fig. 1. 
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FIG. 4. Schematic illustration of the C12 'in' and 'out' collision processes. The dashed line 
represents a condensate atom and the solid lines represent thermal atoms. 
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FIG. 5. The condensate [Nc], noncondensate (iV) and total (N) number of atoms as a function 
of time following a momentum quench. The trap has a frequency of vq = 187 Hz, and contains 
5 X 10^ ^^Rb atoms at an initial temperature of Tq = 200 nK. 




FIG. 6. The rms radius, r = \/{r^, as a function of time following a momentum quench; r„ 
is for the thermal cloud and fc for the condensate. 
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I, and the displacement between the mean condensate 
as a function of time. The initial conditions are given in the 
text. The dashed curves are in the absence of collisions, while the solid curves include both C12 
and C22 collisions. 



FIG. 7. The center of mass position, Zc 
and noncondensate positions, Zc - 
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FIG. 8. Angular displacements of the condensate (a) and thermal cloud (b) as a function of 
time. The points are the results of the simulation and the solid lines are three-mode fits to the 
data. The results are for the trap studied experimentally in Ref. |23|. 
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FIG. 9. As in Fig. | but for T = 185 nK. 
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FIG. 10. As in Fig. | but for T = 296 nK. 
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FIG. 11. Density contour plots of the condensate (dotted lines) and thermal cloud (solid 
lines) for a sequence of equally spaced time intervals spanning approximately one-half period of 
the scissors mode oscillation. The earliest time is at the top of the figure. The dashed horizontal 
line denotes the symmetry axis of the trapping potential and the straight solid line indicates the 
major axis of the condensate at each time step. 
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FIG. 12. Frequency (a) and damping rate (b) of the scissors modes for a variable total number 



of atoms, intended to simulate the experiments in Ref. |23|. The condensate mode is indicated by 
open (theory) and solid (experiment) circles. The open squares in (b) show the calculated average 
damping rate of the two thermal cloud modes, while the solid squares are the corresponding 
experimental values. 
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FIG. 13. The ratio of the moment of inertia of the trapped gas to the rigid body moment of 
inertia as a function of temperature. The solid circles are the result of Eq. ( pO| ) while the open 
circles are obtained using Eq. (|35 
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